home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / astronomy / skyview / !SkyView / c / Ecliptic < prev    next >
Encoding:
Text File  |  1993-08-26  |  4.0 KB  |  125 lines
  1. /********************************************************/
  2. /*                   --  Ecliptic  --                   */
  3. /*          Sundry utility functions concerning         */
  4. /*            ecliptic latitude and longitude.          */
  5. /********************************************************/
  6. #include <math.h>
  7.  
  8. #include "os.h"
  9. #include "coords.h"
  10. #include "menu.h"
  11. #include "sv_header.h"
  12. #include "ecliptic.h"
  13.  
  14. /* Set (arbitrary) declination beyond which object is   */
  15. /* considered to be at the celestial pole.              */
  16. #define DEC_LIMIT 1.57
  17.  
  18. /* Set convergence criteria for solving Kepler's eqn.   */
  19. #define MAX_ITER 10         /* Max number of iterations.*/
  20. #define KEP_CONVERGE 1.0E-6 /* Max change in last iter. */
  21.  
  22. /*------------------------------------------------------*/
  23. /* Function returning the obliquity of the ecliptic at  */
  24. /*  the instant implied by T (No. of Julian centuries   */
  25. /*   from noon GMT on 0/1/1900).  Result in radians.    */
  26. /*------------------------------------------------------*/
  27.   double ecliptic_obliq(double T)
  28. {
  29.   return ( 23.452 - 0.013 * T )*DblCONV;
  30. }
  31.  
  32. /*------------------------------------------------------*/
  33. /*   Function to calculate the elements of the sun's    */
  34. /*     "orbit" at the instant defined by T (Julian      */
  35. /*    centuries from 0/1/1900).  Angles in radians.     */
  36. /*        Coeffs from P. Duffett-Smith p. 116.          */
  37. /*------------------------------------------------------*/
  38.   void ecliptic_sunorbit(double T, double *Lptr, double *Mptr, double *eptr)
  39. {
  40.   double T2 = T * T;
  41.  
  42. /* Mean latitude: */
  43.   *Lptr = 2.7969668E2 + (3.600076892E4 * T) + (3.025E-4 * T2);
  44.   *Lptr = fmod(*Lptr, 360.0)*DblCONV;
  45.  
  46. /* Mean anomaly:  */
  47.   *Mptr = 3.5847583E2 + (3.599904975E4 * T) - (1.54E-4 * T2) -
  48.           (3.3E-6 * T * T2);
  49.   *Mptr = fmod(*Mptr, 360.0)*DblCONV;
  50.  
  51. /* Eccentricity:  */
  52.   *eptr = 1.675104E-2 - (4.18E-5 * T) - (1.26E-7 * T2);
  53.  
  54.   return;
  55. }
  56.  
  57. /*------------------------------------------------------*/
  58. /* Function to solve Kepler's Equation.  Follows method */
  59. /*  given in NAO Technical Note Section 5 (Dec 1978 &   */
  60. /*    Feb 1981) p7.  Returns TRUE if E converged OK.    */
  61. /*------------------------------------------------------*/
  62.   BOOL ecliptic_kepler(double M, double e, double *Eptr)
  63. {
  64.   int iter_count = 0;
  65.   double denom, newE;
  66.   BOOL done;
  67.  
  68.   *Eptr = M + e * sin(M);
  69.   do
  70.   {
  71.     iter_count += 1;
  72.     denom = 1 - e * cos(*Eptr);
  73.     if (denom == 0.0) return FALSE;
  74.     newE = *Eptr + (M + e*sin(*Eptr) - *Eptr)/denom;
  75.     done = fabs(*Eptr - newE) <= KEP_CONVERGE;
  76.   } while (!done && iter_count < MAX_ITER);
  77.  
  78.   return (done);
  79. }
  80.  
  81. /*------------------------------------------------------*/
  82. /* Function to convert ecliptic latitude and longitude  */
  83. /* to right ascension and declination (all in radians). */
  84. /*------------------------------------------------------*/
  85.   void ecliptic_radec(REAL elatit, REAL elongit, double T,
  86.                       REAL *raptr, REAL *decptr)
  87. {
  88.   double beta;      /* Ecliptic latitude.               */
  89.   double lambda;    /* Ecliptic longitude.              */
  90.   double epsilon;   /* Obliquity of the ecliptic.       */
  91.   double alpha;     /* Right ascension.                 */
  92.   double delta;     /* Declination.                     */
  93.   double x, y;
  94.   double arg;
  95.  
  96.   beta   = (double)elatit;
  97.   lambda = (double)elongit;
  98.  
  99.   epsilon = ecliptic_obliq(T);
  100.   arg     = sin(beta)*cos(epsilon) + cos(beta)*sin(epsilon)*sin(lambda);
  101.   if (arg < -1.0) arg = -1.0;
  102.   if (arg >  1.0) arg =  1.0;
  103.   delta   = asin(arg);
  104.  
  105. /* Special case when object is very close to celestial  */
  106. /* poles.                                               */
  107.   if (fabs(delta) > DEC_LIMIT)
  108.   {
  109. /* Set right ascension to an arbitrary value (0.0).     */
  110.     *raptr  = ZERO;
  111.     *decptr = (REAL)delta;
  112.     return;
  113.   }
  114.  
  115.   y = cos(beta)*sin(lambda)*cos(epsilon) - sin(beta)*sin(epsilon);
  116.   x = cos(beta)*cos(lambda);
  117.   alpha = atan2(y,x);
  118.   if (alpha <  0.0      ) alpha += 2.0*DblPI;
  119.   if (alpha >= 2.0*DblPI) alpha  = 0.0;
  120.  
  121.   *raptr  = (REAL)alpha;
  122.   *decptr = (REAL)delta;
  123.   return;
  124. }
  125.